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An adjoint-based methodology for design optimization of unsteady turbulent flows on 
dynamic unstructured grids is described. The implementation relies on an existing unsteady 
three-dimensional unstructured grid solver capable of dynamic mesh simulations and discrete 
adjoint capabilities previously developed for steady flows. The discrete equations for the 
primal and adjoint systems are presented for the backward-difference family of time- 
integration schemes on both static and dynamic grids. The consistency of sensitivity 
derivatives is established via comparisons with complex-variable computations. The current 
work is believed to be the first verified implementation of an adjoint-based optimization 
methodology for the true time-dependent formulation of the Navier-Stokes equations in a 
practical computational code. Large-scale shape optimizations are demonstrated for 
turbulent flows over a tiltrotor geometry and a simulated aeroelastic motion of a fighter jet. 


Nomenclature 


a,b,c,d 

Temporal coefficients 

n 

m x xm x block diagonal rotation matrix 

C 

Aerodynamic coefficient 

n 

3x3 rotation matrix 

C T 

Rotor thrust coefficient 

H GCL 

m q x m q diagonal GCL residual matrix 

D 

Vector of design variables 

s 

Control volume surface area 

E 

Total energy per unit volume 

T 

4x4 transform matrix 

F 

Flux vector 

t 

Time 

F,,F y 

In viscid and viscous flux vectors 

u,v,w 

Cartesian components of velocity 

/ 

Cost function 

V 

Control volume 

/>•* 

General functions 

V 

m q x m q diagonal matrix of cell volumes 

G 

Grid operator 

w 

3x1 face velocity vector 

i 

sR 

X 

m x X 1 vector of grid coordinates 

i, j,k,n 

Indices 

X 

3x1 position vector 

in 

Quantity at initial conditions 

X 

Independent variable 

J 

Number of cost function components 

x,y,z 

Cartesian coordinate directions 

K 

m x Xm x linear elasticity coefficient matrix 

£ 

Perturbation 

L 

Lagrangian function 

A / 

m q X 1 flowfield adjoint variable 

m x 

Size of vector X 

A g 

m x X 1 grid adjoint variable 

m 

Size of vector Q 

0 

Rotor blade collective setting 

N 

Number of time levels 

P 

Density 

n 

Outward-pointing normal vector 

T 

3x1 translation vector 

P 

Pressure, also cost function exponent 

X 

m x X 1 translation vector 

Q 

m q X 1 vector of volume-averaged 

CO 

Cost function component weight 


conserved variables 

'P 

Rotor azimuth 

q 

m q xl vector of conserved variables 

oo 

Quantity at freestream conditions 

R 

m q xl vector of spatial undivided residuals 

* 

Target quantity 
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Introduction 

As computational fluid dynamics (CFD) tools become more efficient, accurate, and robust, their role 
in the analysis and design of new aerospace configurations continues to increase. Computational methods 
have already become a major integrated component of industrial practices. The use of CFD has been 
traditionally confined to the steady regime; however, with recent algorithmic improvements and the 
persistent growth of computational power, CFD methods have begun to make substantial in-roads in 
simulating unsteady flow phenomena. Target applications for these methods are widely abundant; typical 
examples might include the prediction of aeroelastic characteristics, maneuvering flight conditions, six 
degree-of-freedom simulations, specified motion problems, or flow control simulations, among many 
others. 

In recent years, steady-state CFD methods have been targeted for use in automated design optimization 
frameworks. In gradient-based design approaches, one of the major challenges is to obtain sensitivity 
information for the flowfield at a reasonable cost. Conventional black-box finite difference methods 111 
suffer from well-known step-size limitations and incur a computational expense that grows linearly with 
the number of design variables. Forward, or direct, differentiation methods 121 and techniques based on the 
use of complex variables 131 mitigate the step-size limitation but still suffer from excessive cost in the 
presence of many design variables, as is often the case with aerodynamic design applications. 

Adjoint methods provide a powerful alternative for aerodynamic sensitivity analysis. In this approach, 
the sensitivities of an objection function are determined through the solution of an auxiliary, or adjoint, set 
of equations. Adjoint methods may be further categorized into either continuous or discrete approaches, 
depending on the order in which the governing equations are differentiated and discretized. The discrete 
approach allows one to account for mesh variation as well; a second adjoint system can be solved to 
linearize the relationship between the design variables and the mesh operator as described in Ref. [4], The 
principal advantage of the adjoint approach is that the computational cost is independent of the number of 
design variables; a rigorous sensitivity analysis for hundreds of variables can be performed at a cost 
equivalent to the solution of the governing equations themselves. For examples of the use of such methods, 
see the references cited in Ref. [5], 

The role of adjoint-based methodologies in mesh adaptation strategies should also be noted. Whereas 
many traditional mesh adaptation schemes rely on heuristic connections between solution gradient 
information and local mesh spacing requirements, the adjoint equations establish a rigorous mathematical 
connection between solution accuracy and the computational grid. The approach has proven quite powerful 
and has enjoyed success where traditional feature-based approaches have consistently failed. Ref. [6] 
provides a review of recent applications and an extensive list of references on the subject. 

Some recent examples of adjoint-based strategies for unsteady aerospace applications are given in 
Refs. [7]-[ 13], The goal of the current work is to extend the time-dependent adjoint formulation for static 
grids introduced in Ref. [13] and the steady-state discrete adjoint capability developed in Refs. [4] and 
[ 14]-[ 1 8] to the three-dimensional time-dependent Euler and Reynolds-averaged Navier-Stokes equations. 
The present approach and implementation are valid for unsteady flows on various grids including static 
grids, dynamic grids undergoing rigid motion, and general morphing grids governed by a mesh deformation 
scheme based on a linear elasticity analog. The work is believed to be the first verified implementation of 
an adjoint-based optimization methodology for the true time-dependent formulation of the Navier-Stokes 
equations in a practical computational code. In the following sections, the unsteady governing equations are 
presented as well as various mesh motion strategies. This is followed by the derivation of the discrete 
adjoint equations for the flowfield and mesh, including details concerning their implementation. Examples 
demonstrating the discrete consistency of the implementation and applications of the design optimization 
framework to large-scale problems are also shown. 

The Flowfield Equations 
Governing Equations and Solution Method 

Using the approach outlined in Ref. [19], the unsteady Euler and Navier-Stokes equations may be 
written in the following form for both moving and stationary control volumes: 

qrfV + <£ (F f -F v )-fWS = 0, (1) 

at Jv J dV 
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where V is the control volume bounded by the surface dV . The vector q represents the conserved 
variables for mass, momentum, and energy, and the vectors F,- and F v denote the inviscid and viscous 
fluxes, respectively. Note that for a moving control volume, the inviscid flux vector must account for the 
difference in the fluxes due to the movement of control volume faces. Given a flux vector F on a static 
grid, the corresponding flux F,- on a moving grid can be defined as F ; =F — q(W-n), where W is a local 
face velocity and n is an outward-pointing unit face normal. 

By defining a volume-averaged quantity Q within each control volume, 


f qdV 

Q = E , 

V 


(2) 


the conservation equations take the form 


9(QVQ , 

dt 



(F; - F v ) • ndS = 0 , 


(3) 


where the conserved variables and inviscid flux vectors are defined as Q = [p, pu, pv, pw, E\ and 


p{u-w x ) 1 p(v-w y ) r p(w—w z ) 

pu(u-W x ) + p pu(v-Wy) pu{w-W z ) 

F ; = pv(u-W x ) i+ pv(v - W y ) + p j+ pv{ vv - W z ) k. (4) 

pw(u-W x ) pw(v-W y ) pw(w-W z )+ p 

_(E+ p)(u-W x ) + W x p_ ^E + p){v-W y ) + W yP L (E + p)(w-W z ) + W z p_ 


The viscous flux vector F v is not explicitly shown here. The equations are closed with the perfect gas 
equation of state and an appropriate turbulence model for the eddy viscosity. Finally, it is worth noting that 
for the special case of a spatially and temporally constant state vector, e.g., Q = (1, 0, 0, 0, 0) , the 
conservation equations reduce to the Geometric Conservation Law (GCL) 1 

^ = (£> Wf idS. (5) 

dt J dV 


In computational practice, the discrete GCL residual is added to the flow equations to preserve a constant 
solution on dynamic grids/ 191 

References [14], [19], [21], and [22] describe the flow solver used in the current work. The code can be 
used to perform aerodynamic simulations across the speed range and an extensive list of options and 
solution algorithms is available for spatial and temporal discretizations on general static or dynamic mixed- 
element unstructured meshes which may or may not contain overset grid topologies. 

In the current study, the spatial discretization uses a finite-volume approach, in which the dependent 
variables are stored at the vertices of tetrahedral meshes. Inviscid fluxes at cell interfaces are computed 
using the upwind scheme of Roe, 1231 and viscous fluxes are formed using an approach equivalent to a finite- 
element Galerkin procedure. For dynamic mesh cases, the mesh velocity terms are evaluated using 
backward differences consistent with the discrete time derivative; this makes the spatial and GCL residuals 
dependent on grids at previous time levels. The eddy viscosity is modeled using the one-equation approach 
of Spalart and Allmaras. 1241 Massively parallel scalability is achieved through domain decomposition and 
message passing communication. 

An approximate solution of the linear system of equations formed within each time step is obtained 
through several iterations of a multicolor Gauss-Seidel point-iterative scheme. The turbulence model is 
integrated all the way to the wall without the use of wall functions and is solved separately from the mean 
flow equations at each time step with a time integration and linear system solution scheme identical to that 
employed for the mean flow equations. 


The Grid Equations 

The general grid equations can be defined in the form G” (X,D) = 0 , where X is the mesh (meshes at 
several time levels may be involved), D is the vector of design variables, and n denotes the time level and 
indicates that the grid operator may vary in time. The specific formulations for different grid motions are 
introduced below. 
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Grids Undergoing Rigid Motion 

For problems in which rigid mesh motion is required, the motion is generated by a 4x4 transform 
matrix, T , as outlined in Ref. [19]. This transform matrix enables general translations and rotations of the 
grid according to the relation 

x = Tx° , (6) 

o o o o r •• t 

which moves a point from an initial position x =(x ,y ,z ) to its new position x = (x, y, z.) : 

* ^11 R \2 ^13 T x x ° 

y _ R 2l R 22 R 23 T y V° 

z ^31 R 32 R 33 t z z° 

_ij [ o ° ° 1 J i 

In an expanded form, x = IZx 0 + t . Here, the 3x3 matrix TZ defines a general rotation and the vector T 
specifies a translation. The matrix T is generally time-dependent. One useful feature of this approach is 
that multiple transformations telescope via matrix multiplication. This formulation is particularly attractive 
for composite parent-child body motion, where motion of one body is often specified relative to another. 
The reader is referred to the discussion in Ref. [19] for more details. For this formulation, the grid operator 
at time level n is defined as 

G"(X",X°,D) = TZ n X° +x n -X” , (8) 

where X° and X" are the grid vectors at the initial and n th time levels, respectively; TZ" is an m x X m x 
block-diagonal matrix with 3x3 blocks representing rotation and m x being the size of vector X” ; and %" 
is an m x -size translation vector. The matrix 7 Z n and vector x” may explicitly depend on D . 

Deforming Grids 

The simplest example of a deforming grid simulation is a static grid undergoing deformations as a 
result of a shape optimization process. In this case, the grid is not time-dependent and is modeled as an 
elastic medium which obeys the elasticity relations of solid mechanics. An auxiliary system of linear partial 
differential equations (PDE’s) is solved to determine the mesh coordinates after each shape update. 
Discretization of these PDE’s yields a system of equations 

KX = X sulf , (9) 

where 7C represents the elasticity coefficient matrix, X is the vector of grid coordinates being solved for, 
and X sur j- is the vector of updated surface coordinates, complemented by zeros for all interior coordinates. 

The coefficients of the matrix 7C depend on the coordinates of the grid. In the approach followed here, 
the elasticity equations are discretized on the grid corresponding to the initial time level. Thus, the grid at 
the initial level satisfies the nonlinear equations 

K° (X° ,D)X° =X° sutf . (10) 

The material properties of the system are chosen based on the local cell geometry and proximity to the 
surface, and the system is solved using a preconditioned GMRES algorithm. For further details on the 
approach, see Refs. [16], [19], and [25]. 

For static grid formulations, the only grid operator used at all times is 

G(X,D) = X sur j -7CX , (11) 

where X sur f may explicitly depend on D . There are situations in which time-dependent deforming grids 
are required, including aeroelastic deflections of the surface, for which the rigid motion as described in the 
previous section is not valid. Instead, a morphing mesh formulation is used. In this approach, the linear 
elasticity equations given by Eq. (9) are solved at each time level with the matrix 7C = 7C computed at the 
initial time level and fixed throughout the time evolution; the vector X” represents the current body 
positions. For morphing grids, the operator at time level n is defined as 

G" (X" , D) = X n surf - 7C°X" . (12) 

When the surface motion is governed by the rigid motion relations given by Eq. (6), X" () y can be 
further specified as X" u ,y = 1Z" xP sw j + x" . 
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Cost Functions 

The steady-state adjoint implementation described in Refs. [4] and [14]-[18] permits multiple objective 
functions and explicit constraints of the following form, each containing a summation of individual 
components: 

fi = i.<0j(Cj-Cj) p ' . (13) 

7=1 

Here, COj represents a user-defined weighting factor; C ; is an aerodynamic coefficient such as the total 
drag or the pressure or viscous contributions to such quantities, the (*) superscript indicates a user-defined 
target value of C j ; and p ; is a user-defined exponent chosen so that f] is a convex functional. The user 
may specify computational boundaries to which each component function applies. The index i indicates a 
possibility to introduce several different cost functions or constraints, which may be useful if the user 
desires separate sensitivities, e.g., for lift, drag, pitching moment, etc. 

For the unsteady formulation, similar general cost functions /)•” are defined at each time level n . The 
integrated cost function f t is defined as a discrete time integral over a certain time interval [ tj , tf ]: 

N? 

f i = Y j f" A/, (14) 

n=N\ 

12 12 * 

where time levels Nj and Nj correspond to tj and t t , respectively. The user now supplies time intervals 
over which the cost functions are to be used. 


Derivation of the Time-Dependent Adjoint Equations 

To derive the time-dependent form of the adjoint equations, the methodology developed in Ref. [13] is 
used. The governing equations given by Eq. (3) are rewritten as 

-^2 + R = 0, R = ^)(F ; -F v )-n(/5. (15) 

dV 

Using a first-order backward difference (BDF1) in time, the equations can be evaluated at time level n as 
follows: 

V” Q ~ Q — + R" + 'TZq CL Q "~ 1 = 0 . (16) 

At 

Here, V” and are m q x m q diagonal matrices, m q is the length of vector Q" , the GCL is 

discretized in a consistent fashion as 


± ( v' ! -v"- 1 ) = ^ cl . 

At 


(17) 


and R” is the spatial undivided residual. Recall that R" and 'R-qcl depend on grids at the current and 
previous time levels. Note also that although the BDF1 scheme has been shown here for the sake of 
simplicity, the derivations for higher-order temporal schemes are similar and included as an appendix. 

The discrete adjoint-based optimization methodology is based on the method of Lagrange multipliers, 
which is used to enforce the governing equations as constraints. For the sake of simplicity in the following 
presentation, a single cost function is assumed and therefore the index i is omitted. For the time-dependent 
equations, the Lagrangian functional is defined as follows: 


N 


N 


L(D.Q. X. A f , A, ) = X / " 4 ' + Z[A/’ 


V 


n = 1 


72=1 


Q'-Q 

At 


72—1 


-+R” +Hq CL Q 


72—1 


At 


( 18 ) 


N 

£[' 


72=1 


k] 

7 G"A/ + 

+ 

o 

V 

A °/’ 

r R'" |a/ + 




V 


J 



G°At 


where f n = 0 for n < N l and n > N 2 ; G" = 0 are the grid equations at time level n ; A" and A" are 
vectors of Lagrange multipliers associated with the flow and grid equations at time level n , respectively; 
D is a vector of design variables; and R m = 0 is the initial condition for the flow equations. 
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The Lagrangian is differentiated with respect to D , assuming that /" depends on Q" , X" , and D ; 
R'" depends on Q° , X° , and D ; R" depends on Q" , X" , X" -1 , and D ; and 7 Z£ CL depends on X" , 
X" -1 , and D. Regrouping terms to isolate the coefficients of 3Q” /dD and equating the coefficients to 
zero yields the final form of the adjoint equations for the flowfield: 


— (v 

At ' 


— |V"A'}-V 


71+1 a 71+1 ' 


l / 


3R" 

3Q" 


IT 


\ n A n+1 - _ 

A/ + /C GC L a / - 


jr 

3Q" 


1 T 


, for 1 < n < N 


~—V l A l f + 
At J 


3R" 


-\T 


3Q 


o 


A/ +7^-gclA/ -- 


r 

3Q° 


it 


, for the initial time level. 


(19) 


(20) 


where A /V 4 ’ 1 = 0 . Collecting the coefficients of dX"/rJD and equating them to zero leads to similar adjoint 
equations for the grid. Assuming that the grid operator at time level n , G" , depends on X" , X° , and D , 
the grid adjoint equations are defined as 

lT 


1 

QJ 

o 

a 

i 

T 

A" - 

r 1 

T 

3V Q" -Q" -1 

ax" _ 


_ax"_ 

1 

ax" At 


\ n 

A / 


( 21 ) 
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3R 


3G 

ax° 


-T 


A ( > = 
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z 

71=1 L 


dG n 

ax° 


k=0 L 

. lT 


ax" 


n+k -VT7 n+k 

! d/< -GCL Q«+t:-l 


-T 


ax" 


A n f +k , for 1 < n < N ; 
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a/ 0 

ax° 
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aR" : 
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~\T 


A / + 


aR | diz (;CL q 0 


~\T 


ax u 


ax u 
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(22) 


The specific form of these equations will be discussed in subsequent sections. With the adjoint coefficients 
satisfying the flowfield and grid adjoint equations, the sensitivity derivatives are calculated as follows: 


dL 

dD 


N 




aR" | dTV GCL q„_i 


3D 


3D 


0 


3D 

t 3R'" - 

3D 


+ 

kT SG "l 


L a D 

/ 


At 


(23) 


T 3G° 


At . 


Implementation 


Flowfield Adjoint Equations 

The implementation and solution of Eqs. (19) and (20) are based largely on the steady-state strategies 
described in Refs. (4] and [ 14]-[ 1 8] . In this manner, a great deal of software development effort is avoided 
since the steady and unsteady equations share many similar terms, namely the details of the spatial 
discretization. However, some fundamental differences in the implementation must be addressed for time- 
dependent problems. 

Implications of Reverse Time Integration 

While the discrete solution Q" for Eq. (3) is determined by marching forward in physical time from 
n = 0 to n = N , due to the nature of the adjoint equations and their boundary conditions, the solution for 
A" must instead be initiated from n = N and proceed backwards in physical time. Since Eqs. (19) and 
(20) involve the linearizations 3R"/3Q and df"/dQ, the flow solution Q" at all time levels must be 
available during the reverse integration. 

In practice, the most straightforward approach to meeting this requirement is to store Q" to disk for all 
n during the solution of Eq. (16). In this case, the storage cost is significant, but the primary advantage is 
ease of implementation. This is the approach used for the current study. For problems where the mesh is 
changing in time, the grid point coordinates and associated speeds are also stored. While these mesh-related 
values could be recovered by performing the mesh motion in reverse, ease of the full storage 
implementation has been favored. 
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Solution Strategy 

As described in Ref. [19], each solution vector Q” is determined through a dual time-stepping 
procedure. In this approach, a sequence of subiterations is performed within each physical time step. The 
procedure relies on an approximate linearization of the discrete residual combined with a pseudo-time term 
to achieve a scheme directly analogous to that used in Ref. [21] for steady flows. The same subiterative 
strategy is employed for the time-dependent adjoint equations, following an approach similar to that 
outlined in Ref. [17]. The Jacobian matrix used to relax the adjoint system is constructed once at each time 
step n based on the value of Q” , and does not change during the subiterative procedure. 

A requirement for performing adjoint solutions is that the iteration scheme be linearly stable. It has 
been observed in some cases, more often for unsteady problems than for steady ones, that linear stability is 
not satisfactory. Suggested explanations [181[26I,[271 vary from physical instabilities to instabilities of the 
numerical schemes involved. The generalized conjugate residual (GCR) scheme described in Ref. [28] has 
been used to wrap the multicolor Gauss-Seidel iteration as well as the temporal subiterative procedure. This 
approach has been found to work well in stabilizing otherwise problematic iterations. 

Data Storage 

For three-dimensional dynamic grid simulations using a one-equation turbulence model, the reverse 
time integration and solution techniques outlined above require the storage of twelve floating-point 
variables per grid point at each time step: six flowfield variables, three mesh coordinates, and three mesh 
velocities. For large-scale problems involving many time steps, this strategy can easily result in a storage 
requirement on the order of terabytes of disk space. Strategies for circumventing storage limitations have 
been suggested in the literature; [9I ’ [29I,[301 these may be the focus of future investigations once an initial 
capability has been established. 

In the current implementation, each processor is responsible for reading and writing its local solution 
for the entire time history to a unique file on disk. Since each file may contain several gigabytes of data, 
requiring several hundred processors to parse sequential-access files at each time step can be very 
inefficient. For this reason, direct-access files are used so that the file pointer can be immediately placed at 
the record of interest. It has been found that this approach can decrease the time required for disk I/O by as 
much as two orders of magnitude on large cases. The use of asynchronous file I/O was also examined, 
although it is not currently being used. 

Grid Adjoint and Sensitivity Equations 

Depending on the nature of the grid operator G and the design variables D , the grid adjoint and 
sensitivity equations may need to be solved at each time level n , once at n = 0 , or not at all. If solutions at 
each time step are required, they are performed at the completion of each step of the adjoint solver, rather 
than subsequently performing additional loops over the entire range of time levels. In this manner, Q" is 
the only solution vector which must be stored for all n , while A" and A” may be discarded when no 
longer needed. 

The predominant challenge in the discretization and solution of Eqs. (21)-(23) is the infrastructure 
required to simultaneously manage data from several time levels. Inspection of Eqns. (42)-(44) in the 
appendix that are higher-order analogs to Eq. (21) shows that for a given time step n , the solution for A" 
may depend on values of Q from adjacent time levels both prior and subsequent to level n . Values of A f 
must also be available at time level n as well as later time levels. Moreover, this complexity increases with 
the temporal order of the scheme. 

The summation term in Eq. (21) is ultimately due to the dependency of the mesh speeds on grid 
coordinates at multiple time levels, according to the BDF scheme being used. Rather than linearizing R 
and V, GCL at several time levels with respect to the grid coordinates at the current time level as indicated 
in the summation, an inverse approach more amenable to the existing implementation of the spatial 
linearizations is used. The residual at time level n is linearized with respect to the grid coordinates at every 
time level in the temporal stencil by seeding the linearizations with the appropriate BDF coefficient. The 
results are then stored temporarily for use in evaluating the summation term at subsequent time levels 
within the stencil, after which the linearizations are discarded. 
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Verification of Adjoint Implementation 

To verify the accuracy of the implementation, comparisons are made with results generated through an 
independent approach based on the use of complex variables. This approach was originally suggested in 
Refs. [31] and [32], and was first applied to a Navier-Stokes solver in Ref. [3], Using this formulation, an 
expression for the derivative of a real-valued function fix) may be found by expanding the function in a 
complex-valued Taylor series, using an imaginary perturbation ie : 

3/ Im{f(x + ie)\ 2 .... 

— = — — - + 0(e ). (24) 

dx £ 

The primary advantage of this method is that true second-order accuracy may be obtained by selecting step 
sizes without concern for subtractive cancellation error typically present in real-valued divided differences. 
Through the use of an automated scripting procedure outlined in Ref. [33], this capability can be 
immediately recovered at any time for the baseline flow solver. For computations using this method, the 

in 

imaginary step size has been chosen to be 10 . For each verification test, all equations sets are converged 

to machine precision for both the complex-variable and adjoint approaches. When used, the elasticity 
matrix K is assumed to be constant throughout the verification. 


Static Grid 


Test Case 

The first test case is used to verify the implementation for unsteady flows on static grids. For this 
example, fully turbulent flow over the ONERA M6 wing [ 41 shown in Fig. 1 is considered. The grid 
contains 16,391 nodes and 90,892 tetrahedral elements and 16 processors are used for the simulation. The 
freestream Mach number is 0.3, the angle of attack is 1°, and the Reynolds number is 1 million based on the 
mean aerodynamic chord (MAC). The simulation is initiated from freestream conditions Q°° , which leads 
to R'" = Q°° - Q° . The solution is advanced 5 physical time steps using a nondimensional At of 0.1. 
While this coarse spatial resolution, relatively large time step, and brief duration of the simulation are not 
sufficient to resolve the flow physics of the problem, they are adequate to evaluate the discrete consistency 
of the implementation. 

Design Variables 

For this test, two general classes of design variables are used. The first class of variables is comprised 
of global parameters unrelated to the computational grid. These variables include parameters such as the 
freestream Mach number and angle of attack. Such variables are useful in verifying the implementation of 
the flowfield adjoint equation, as the terms in Eq. (23) associated with these parameters are generally trivial 
to implement or identically zero, and solution of the mesh adjoint equations is not required. 

The second class of design variables provides general shape control of the configuration. The 
implementation allows the user to employ a geometric parameterization scheme of choice, provided the 
associated surface grid linearizations are available. For all examples in the current study, the grid 
parameterization approach described in Ref. [35] is used. This approach can be used to define general 
shape parameterizations of existing grids using a set of aircraft-centric design variables such as camber, 
thickness, shear, twist, and planform parameters at various locations on the geometry. The user also has the 
freedom to associate two or more design variables to define more general parameters. In the event that 
multiple bodies of the same shape are to be designed, the implementation allows for a single set of design 
variables to be used to simultaneously define such bodies. In this fashion, the shape of each body is 
constrained to be identical throughout the course of the design. 


Grid Adjoint Equation 

For this case, there is only one grid operator, G(X,D) = X w y -XIX , which does not depend on time. 
As a result, the grid adjoint equation can be recast as 
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and the sensitivity derivative is 

N V r 3R" 


3L 
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Computational Results 

The test has been performed using the BDF1 scheme and all other time-integration schemes described 
in the appendix, and results are listed in Table 1. Sensitivity derivatives of the lift coefficient at the final 
time step with respect to the angle of attack and a camber variable located at the midspan of the wing are 
shown. The results for the adjoint implementation exhibit excellent agreement with the complex-variable 
approach, differing at most in the fifteenth digit. 


Rigidly-Moving Grid 

Test Case 

The next test case is used to verify the implementation for rigidly moving meshes. For this case, the 
grid and freestream conditions and computational environment are identical to those described for the 
previous test; however, the mesh is now subjected to an oscillatory pitch-plunge motion based on the rigid 
mesh transform approach outlined earlier. The nondimensional pitching and plunging reduced frequencies 
are 0.5 and 0.1, respectively. The pitching amplitude is 5° and takes place about a vector normal to the 
symmetry plane located 0.47 MAC from the wing root leading edge. The amplitude of the plunging motion 
is 0.38 MAC. The baseline wing position at t = 0 is as shown in Fig. 1. As in the prior test, the simulation 
is initiated from freestream conditions R'" = Q°° - Q° and is advanced 5 physical time steps using a 
nondimensional At of 0. 1 . 


Design Variables 

The design variables for the current test include those described above for the static grid example, as 
well as a third class of parameters governing the rigid motion procedure described above. These include 
translation and rotation frequencies, amplitudes, and directional vectors, as well as centers of rotation. 


Grid Adjoint Equation 

For this test case, the following grid operators are used: at the zeroth time level, the grid is either 
unchanged or governed by the elasticity equations G°(X°,D) = X^y -7C°X° ; grids at other time levels 
are governed by the rigid motion equation G" (X" ,D) = 1Z"X° + z n - X" . 

The grid adjoint equations are given by 

a; lVy[— +^ Q ”HV-f°r '***»■ m 
S Lax"J 3X" A t 1 *-i 3X" 3X" 1 

L J x=(J I- -I 

Under the assumption that the shape does not change ( X° is constant), the sensitivity derivative is given 
by 
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The formulation that would allow shape design is the following: 
nT ( N r } r-vfO 

£[K"] r ."f l» 

V n = 1 


o,3/C „o 


1C’ + 


3X U 


A° = 


A" 

8 


+ 

J 


3X U 


\T r 

3CT 

3X° 


A / + 


3r | dn GCL qc 


-\T 


3X 


and the corresponding sensitivity derivative is 
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Computational Results 

Results for the derivatives of the lift coefficient at the final time step are shown in Table 2 for the 
current case. In addition to the angle of attack and camber variables, derivatives with respect to the rigid 
motion pitching frequency are also shown. The agreement with the complex-variable formulation is 
excellent for each of the time integration schemes considered. 


Morphing Grid 

Test Case 

To evaluate the accuracy of the implementation for morphing grids, the test case used for rigid motion 
described above is repeated with slight modifications. For the current test, the surface grid of the wing is 
moved using rigid motion, while the interior of the mesh is determined using the elasticity relation given by 
Eq. (9). All other input parameters remain unchanged. 

Design Variables 

The current test case uses the same design variables as the rigid motion test case described above. 

Grid Adjoint Equation 

At all time levels, the grids are governed by the elasticity equations G”(X”,D) = X” (/ y -7C°X" , and 
the surface coordinates are governed by the rigid motion equation X” u) y = TV' X ( s * u ,y + z" . 

The grid adjoint equations are given by 



The sensitivity derivative is 



Two observations can be made. First, note that in the absence of any surface motion, i.e., 7 Z" is the 
identity matrix and t" = 0 , the morphing grid formulation is equivalent to the static grid formulation. Also, 
with a constant transformation matrix T applied to all computational boundaries, the morphing and rigidly 
moving grid formulations are equivalent. 

Computational Results 

The results for the current test case are shown in Table 3. Derivatives of the lift coefficient at the final 
time step with respect to each of the design variables exhibit excellent agreement for the adjoint 
implementation and complex variable formulation. 


Large-Scale Design Cases 

Two large-scale design optimization examples are presented. Although the grid motion in both cases is 
prescribed, a more realistic treatment would involve the use of additional coupled computational models 
such as six degree-of-freedom or structural simulations. While such capabilities are available for use with 
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the flow solver, [191 their effects have not been accounted for in the derivation and implementation of the 
adjoint equations. This important development is relegated to future work. 

Both of the example cases shown below have been performed using 128 dual-socket quad-core nodes 
with 3.0 GHz Intel Xeon processors in a fully-dense fashion for a total of 1,024 computational cores. This 
environment has been chosen to maximize computational efficiency for the chosen test problems; 
numerical experiments have shown that the solvers used in the current study scale well in this range for the 
grid sizes selected. 

The computational grid sizes and time steps for the examples presented here have been chosen merely 
to demonstrate optimization capability for typical problems using immediately available resources. Spatial 
and/or temporal refinement could be readily performed if desired. Although the formulation places no 
restrictions on initial conditions, all solutions are started from freestream conditions. The grids have been 
generated using the method of Ref. [36] and the optimizations have been performed using the package 
described in Ref. [37]. 

Tiltrotor Configuration 

The first large-scale example is a three-bladed tiltrotor configuration similar to that used by the V-22 
aircraft and is based on the Tilt Rotor Aeroacoustics Model (TRAM) geometry described in Refs. [38] and 
[39]. The grid used for this computation is designed for a blade collective setting of © = 14° and consists 
of 5,048,727 nodes and 29,802,252 tetrahedral elements. The rotational speed of the rotor is held constant 
at a value corresponding to a tip Mach number of 0.62 in a hover condition. The Reynolds number is 2.1 
million based on the blade tip chord. The physical time step is chosen to correspond to T of rotor azimuth, 
for a total of 360 time steps per revolution. The BDF2 opt formulation outlined in Ref. [40] is used with 10 
subiterations per time step. 

For this test, the prescribed rigid mesh motion consists of 4 initial revolutions of the geometry 
designed to reach a quasi-steady hover condition, followed by five additional revolutions during which a 
90° constant-rate pitch-up maneuver into a forward-flight mode is performed. A more realistic pitch-up 
scenario might consist of many more revolutions; however, the prescribed motion was chosen to keep the 
cost of the computation affordable on the current resources. During the pitch-up phase of the motion, an 
assumed forward-flight velocity profile based on a simple sine function is imposed through the mesh speed 
terms. The schedule for the shaft angle and forward-flight velocity is shown in Fig. 2, where the shaft angle 
is defined to be 0° in the hover condition and 90° in forward flight. The resulting motion is shown in Fig. 
3, where a snapshot of the rotor is shown every 360° during the course of the motion. An isosurface of the 
second invariant of the velocity-gradient tensor, also known as the Q criterion from Ref. [41], at the time 
step corresponding to T = 1440° is shown in Fig. 4. The tip vortex system is maintained for 2 to 3 
revolutions of the rotor. 

The objective function for the current test case is to maximize the rotor thrust coefficient over the time 
interval corresponding to the pitch-up maneuver, 1441° < T 1 < 3240° : 

3240 

/= ^ (C£-0.1) 2 A/. (34) 

n=1441 

Here, the target thrust coefficient value of 0.1 has been chosen to sufficiently exceed the baseline thrust 
profile shown as the solid line in Fig. 5. After the first four rotor revolutions, the thrust coefficient has 
reached a quasi-steady value of approximately 0.015, which is in good agreement with experimental data 
given in Refs. [38] and [39]. The thrust coefficient shows a discontinuous behavior at the impulsive start of 
the pitch-up motion (n = 1441) and gradually decreases to a lower constant value in the forward-flight 
condition. A subtle 3/rev oscillation in the thrust coefficient during the pitch-up maneuver can also be seen. 

The surface grid has been parameterized as described in Ref. [42]. This approach yields a set of 44 
active design variables describing the thickness and camber of the blade geometry as shown in Fig. 6; 
thinning of the blade is not allowed. In addition, a single twist variable is used to modify the blade 
collective setting during the design. 

The convergence history for six design cycles is shown in Fig. 7. The optimizer quickly reduces the 
value of the objective function over the first two design cycles, after which further improvements are 
minimal. Closer inspection of the design variables indicates that the majority of values have reached their 
bound constraints, preventing any further reduction in the objective function. The final thrust coefficient 
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profile is included as the dashed line in Fig. 5. Cross-sections of the baseline blade geometry are compared 
with the optimized geometry in Fig. 8. The optimization has increased the camber of the blade across the 
span, as well as the blade collective setting. 

The cost of each solution to the unsteady flow and adjoint equations for the current example is 
approximately 3.5 and 10.5 wallclock hours, respectively; however, due to frequent file I/O, this estimate 
varies with file system load. The optimization procedure requires twelve calls to the flow solver and six 
calls to the adjoint solver, for a total runtime of approximately 4.5 days of wallclock time or 1 10,000 CPU 
hours. The disk storage required for one complete flow solution is approximately 1.5 terabytes. 

Fighter Jet with Simulated Aeroelastic Effects 

The second example uses a deforming grid approach to simulate aeroelastic motion of the modified F- 
15 fighter jet configuration described in Ref. [43] and known as NASA research aircraft #837, shown in 
Fig. 9. The computational model assumes half-plane symmetry in the spanwise direction. The grid consists 
of 4,715,852 nodes and 27,344,343 tetrahedral elements and includes detailed features of the external 
airframe as well as the internal ducting upstream of the engine fan face and plenum/nozzle combination 
downstream of the turbine. For the current test, the freestream Mach number is 0.90, the angle of attack is 
0° , and the Reynolds number based on the MAC is 1 million. The static pressure ratio at the engine fan 
face is set to 0.9 and the total pressure ratio at the plenum face is ramped linearly from 1.0 to its final value 
of 5.0 over the first 50 time steps. 

The prescribed grid motion consists of 5 Hz 0.3° oscillatory rotations of the canard, wing, and tail 
surfaces about their root chordlines, with the wing oscillations 180° out of phase with the canard and tail 
motion. In addition, the main wing is also subjected to a 5 Hz oscillatory twisting motion whose amplitude 
decays linearly from 0.5° at the wing tip to 0° at the wing root and takes place about the quarter-chord 
line. This composite motion results in a maximum wing tip deflection of approximately 1.3% MAC as 
shown in Fig. 10. The BDF2 opt scheme is used with 10 subiterations and a physical time step corresponding 
to 100 steps per cycle of grid motion. 

The unsteady lift-to-drag ratio ( L/D ) for the baseline configuration undergoing the specified motion 
for 300 time steps is shown as the solid line in Fig. 11. The L/D behavior begins to exhibit a periodic 
response after approximately 100 time steps. The high-frequency oscillations in the profile are believed to 
be due to a small unsteadiness in the engine plume shown in Fig. 12; this behavior is also present when the 
mesh is held fixed. 

The objective function for the current test case is to maximize LlD on the interval 201 <n< 300 : 

300 2 

/= ^ [(T/T>)' ! -5.°]“a/, (35) 

n=201 

where the target L/D value of 5.0 has been chosen to provide sufficient room for optimization over the 
baseline profile. The surface grids for the canard, wing, and tail have been parameterized as shown in Fig. 
13, resulting in a set of 98 active design variables describing the thickness and camber of each surface. 
Thinning of the geometry is not permitted. 

Convergence of the objective function is shown in Fig. 14. A large reduction in the function is 
obtained after a single design cycle, after which further improvements are minimal due to many of the 
design variables having reached their bound constraints. The final L/D profile is included as the dashed 
line in Fig. 11. The resulting shape changes at various spanwise stations on the canard, wing, and tail are 
shown in Fig. 15, where the vertical scale has been exaggerated for clarity. The design procedure has 
increased the thickness of the wing and canard, as well as the camber across all three elements. Closer 
inspection shows that the trailing edges of each surface have also been deflected in a downward fashion. 

The wallclock times required for single flow and adjoint solutions for the current problem are 
approximately 1 and 1.5 hours, respectively. For the five design cycles shown in Fig. 14, the optimizer 
requires ten flow solutions and five adjoint solutions, or a total wallclock time of approximately 18 hours or 
18,400 CPU hours. The disk space necessary to store a single unsteady flow solution is 136 gigabytes. 
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Summary and Future Work 

A discrete adjoint-based methodology for optimization of unsteady flows governed by the three- 
dimensional Reynolds-averaged Navier-Stokes equations on dynamic unstructured grids has been 
formulated and implemented. The methodology accounts for mesh motion based on both rigid movement 
as well as deforming grids. The accuracy of the implementation has been verified using comparisons with 
an independent approach based on the use of complex variables. The methodology has been successfully 
used in a massively parallel environment to perform two large-scale design optimization examples: one for 
a tiltrotor in a pitch-up maneuver into a forward-flight regime and another for a fighter jet with simulated 
aeroelastic effects. 

Although the approach outlined in the current study represents significant progress towards the goal of 
performing routine optimization of unsteady turbulent flows, a number of research areas remain to be 
explored. The extension of the present formulation to overset grid topologies is ongoing and will allow the 
treatment of multiple bodies undergoing large relative motion. Methods aimed at reducing the storage costs 
associated with the flow solution have the potential to drastically reduce disk requirements. Techniques 
based on variable or adaptive time steps as well as alternate time integration schemes should be examined. 
The effects of related computational disciplines such as six degree-of-freedom and structural models should 
also be properly accounted for. Finally, the use of the unsteady flowfield adjoint solution holds tremendous 
potential for performing mathematically rigorous mesh adaptation to specified error bounds. 
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Appendix: Adjoint Equations for Higher Order BDF Schemes 


The high (up to third) order BDF discretizations for the time derivative of a function s are defined as 


ds 

dt 


1 T n , i n —1 . n —2 . , n — 3 

— as +bs -t-cs +ds 
At L 


(36) 


where n is a time level, and the coefficients are given in Table 4. The coefficients listed for the BDF2 opt 
scheme are a linear combination of the BDF2 and BDF3 coefficients taken from Ref. [40]. The resulting 
scheme is second-order accurate but has a leading truncation error term less than that of the BDF2 scheme. 
While usually found to be stable in practice, stability of the BDF2 opt and third-order BDF3 scheme are not 
guaranteed. Discrete conservation laws are defined as 
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Since the morphing grid formulation includes static meshes and rigid motion as special cases, the 
derivation is provided only for this formulation. Taking into account that R" and R q CL are dependent on 
X” -2 and X ,!_3 , the procedure applied to the BDF1 scheme may also be used to derive the following 
adjoint equations for the flowfield: 
for 3 < n < N : 
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The corresponding mesh adjoint equations are obtained as follows. Assuming 1 1 = R A = R w+3 = 0 
and R^cl = r gcl = r gcl = 0 : 
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(44) 


and for initial conditions, R'” = Q°° - Q° : 
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The sensitivity derivative for the higher-order BDF schemes is evaluated using Eq. (23). 


(45) 
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Table 1. Results for static grid test case. “A” denotes adjoint result, “C” denotes complex- variable result. 


Design Variable 

BDF1 

BDF2 

BDF3 

BDF2 opt 

Angle of 
Attack 

A: 0.004249541855867 
C: 0.004249541855867 

A: 0.003734353591935 
C: 0.003734353591935 

A: 0.003687377975335 
C: 0.003687377975335 

A: 0.003708754474661 
C: 0.003708754474661 

Camber 

A: 0.010713047647152 
C: 0.010713047647155 

A: 0.013701437304586 
C: 0.013701437304586 

A: 0.014574974114575 
C: 0.0145749741 14577 

A: 0.014145698047604 
C: 0.014145698047602 

Table 2. Results for rigidly moving grid. 

“A” denotes adjoint result, “C” denotes complex-variable result. 

Design Variable 

BDF1 

BDF2 

BDF3 

BDF2 opt 

Angle of 
Attack 

A: 0.004713138571667 
C: 0.004713138571667 

A: 0.004293218571759 
C: 0.004293218571759 

A: 0.004245785984455 
C: 0.004245785984455 

A: 0.004267302756747 
C: 0.004267302756681 

Pitching 

Frequency 

A: -0.403740396501207 
C: -0.403740396501207 

A: -0.527819225717431 
C: -0.527819225717432 

A: -0.529833595955533 
C: -0.529833595955533 

A: -0.528894917963836 
C: -0.528894917963837 

Camber 

A: 0.011630821689945 
C: 0.01 1630821689944 

A: 0.013925365539211 
C: 0.013925365539206 

A: 0.014291228334440 
C: 0.014291228334428 

A: 0.014071544549783 
C: 0.014071544549783 


Table 3. Results for morphing grid. “A” denotes adjoint result, “C” denotes complex-variable result. 


Design Variable 

BDF1 

BDF2 

BDF3 

BDF2 opt 

Angle of 
Attack 

A: 0.004713528355526 
C: 0.004713528355526 

A: 0.004298221887378 
C: 0.004298221887378 

A: 0.004250753632738 
C: 0.004250753632738 

A: 0.004272205860974 
C: 0.004272205860974 

Pitching 

Frequency 

A: -0.403961428430834 
C: -0.403961428430834 

A: -0.528263525075847 
C: -0.528263525075847 

A: -0.530205775809711 
C: -0.530205775809710 

A: -0.529295291075346 
C: -0.529295291075346 

Camber 

A: 0.011680362720549 
C: 0.01 1680362720548 

A: 0.013922237526691 
C: 0.013922237526686 

A: 0.014268675858452 
C: 0.014268675858435 

A: 0.014055458873064 
C: 0.014055458873058 


Table 4. Coefficients for higher-order BDF schemes. 


Scheme 

a 

b 

c 

d 

BDF2 

3/2 

-2 

1/2 

0 

BDF3 

11/6 

-3 

3/2 

-1/3 

BDF2 opt 

5.08/3 

-2.58 

1.08 

-0.58/3 
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Figure 2. Forward Mach number and shaft angle schedule for TRAM rotor simulation. 



Figure 3. View of TRAM rotor motion. 
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Figure 4. Isosurface of Q criterion for TRAM rotor at = 1440° . 



Figure 5. Thrust for TRAM rotor before and after design optimization. 
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Figure 6. Spanwise blade and design variable locations for TRAM rotor. 
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Figure 7. Objective function history for TRAM rotor. 



Figure 8. Span wise blade cross-sections before and after optimization of TRAM rotor. 



Figure 9. Modified F-15 with engine duct geometry. 
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Figure 12. Cross-section of engine plume contours for modified F-15. 
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Figure 13. Spanwise and design variable locations for modified F-15. 



Figure 14. Objective function history for modified F-15. 


Wing 



Figure 15. Canard, wing, and tail cross-sections before and after optimization of modified F-15 
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